function x = jacobi(A, b, x0, eps, maxiter)
%	Aceasta functie calculeaza solutia sistemului A * x = b.
%Se considera A = N - P; ( N - P) * x = b; N * x = P * x + b;
%fie N1 inversa matricii N; atunci x = N1 * P * x + N1 * b.
%Pentru constructia matricilor N si P se considera matricile 
%D, L si U. D este o matrice diagonala, ce are pe diagonala 
%principala elementele lui A de pe diagonala principala;
%L este o matrice inferior triunghiulara, ce are sub diagonala
%principala elementele lui A de sub diagonala principala( luate 
%cu semn schimbat);  U este o matrice superior triunghiulara, ce  
%are deasupra  diagonalei principale elementele lui A de deasupra 
%diagonalei principale( luate cu semn schimbat).
%Se observa ca A = D - (L + U). Conform metodei Jacobi, N = D si 
%P = L + U. N fiind o matrice diagonala, se poate calcula usor
%matricea N1( daca toate elemntele de pe diagonala principala sunt
%nenule). 
%	Plecand de la formula: x = N1 * P * x + N1 * b, observam ca solutia
%x se poate calcula recursiv. Astfel, se pleaca de la o aproximare
%initiala x0 si se calculeaza o alta solutie aproximativa, mai aproape de
%cea reala, cat timp numarul de iteratii nu depaseste numarul maxiter
%sau cat timp solutia se incadreaza in toleranta eps impusa, adica cat timp
% || A * x - b || > eps .
%
%	Date de intrare:
%		- A  ->  matricea cu care se inmulteste x in sistemul A * x = b;
%		- b  ->  matricea termenilor liberi;
%		- x0 ->  aproximatia initiala a sistemului;
%		- eps->  toleranta( sau precizia) dorita;
%		- maxiter-> numarul maxim de iteratii.
%
%	Date de iesire:
%		- x  ->  solutia sistemului.

[n n] = size(A);
D = L = U = zeros(n, n);

%mai intai, se afla matricile D, L, U

%matricea D este formata din elementele de pe diagonala principala
for i = 1:n
	D(i, i) = A(i, i);
endfor

%matricea U este formata din elementele de deasupra diagonalei principale ale matricei A(luate cu semn schimbat)  si o in rest
for i = 1:n
	for j = i+1:n
		U(i, j) = -A(i, j);
	endfor
endfor

%matricea L este formata din elementele de sub diagonala principala a matricei A, luate cu semn schimbat  si 0 in rest
for i = 1:n
	for j = 1:i - 1
		L(i, j) = -A(i, j);
	endfor
endfor

%conforma metodei Jacobi, N = D, P = L + U

N = D;
P = L + U;

N1 = inv( N);
G = N1 * P;
c = N1 * b;

iter = 0;
x = G * x0 + c;%se calculeaza un x in afara ciclului 

while( iter <= maxiter && ( norm( A * x - b ) >= eps) )
	iter = iter + 1;
	x = G * x + c;
	x0 = x;%in x0 se retine acel x calculat la pasul anterior
endwhile

endfunction